Quantum Monte-Carlo method applied to Non-Markovian barrier transmission 
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In nuclear fusion and fission, fluctuation and dissipation arise due to the coupling of collective 
degrees of freedom with internal excitations. Close to the barrier, both quantum, statistical and non- 
Markovian effects are expected to be important. In this work, a new approach based on quantum 
Monte-Carlo addressing this problem is presented. The exact dynamics of a system coupled to 
an environment is replaced by a set of stochastic evolutions of the system density. The quantum 
Monte-Carlo method is applied to systems with quadratic potentials. In all range of temperature and 
coupling, the stochastic method matches the exact evolution showing that non-Markovian effects 
can be simulated accurately. A comparison with other theories like Nakajima-Zwanzig or Time- 
ConvolutionLess ones shows that only the latter can be competitive if the expansion in terms of 
coupling constant is made at least to fourth order. A systematic study of the inverted parabola 
case is made at different temperatures and coupling constants. The asymptotic passing probability 
is estimated in different approaches including the Markovian limit. Large differences with the exact 
result are seen in the latter case or when only second order in the coupling strength is considered 
as it is generally assumed in nuclear transport models. On opposite, if fourth order in the coupling 
or quantum Monte-Carlo method is used, a perfect agreement is obtained. 

PACS numbers: 24.60-k,25.70. Jj,05.60.Gg 
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I. INTRODUCTION 

To understand nuclear reactions, the dynamics of nu- 
clei is often replaced by few selected collective degrees of 
freedom expected to contain important information on 
the dynamic. This is for instance the case in fusion reac- 
tions where the relative distance and/or mass asymmetry 
is retained [H, 0] . Another example is provided by the fis- 
sion process which is often treated as a trajectory in an 
energy landscape function on different deformation pa- 
rameters [1, H[ . Although the evolution is projected onto 
few variables, other internal degrees of freedom may play 
an important role to understand the onset of dissipa- 
tion or fluctuation phenomena @. To treat these effects, 
the relevant degrees of freedom should be regarded as an 
Open Quantum System coupled to an environment which 
simulates the internal dynamics. 

To include dissipation in collective space, two impor- 
tant simplifications are often made. First, most of cur- 
rent models treating fusion/fission neglect quantum ef- 
fects and consider a classical treatment 6-9]. Such an 
approximation is however expected to be valid only if 
the internal excitation is high and therefore is not ex- 
pected to hold close to or below the Coulomb barrier. As 
it is discussed in ref. (Iol |. a proper treatment of quan- 
tum and dccoherence effects might be crucial in this re- 
gion. Second, when the time scale associated to collec- 
tive dynamics cannot be dissociated from the one of the 
environment, "non-Markovian" (also called "memory") 
effect should also be properly treated [Tl|, 03 ■ Large ef- 
fort is now devoted to account for both quant al and non- 
Markovian effects in nuclear reactions [llMla] and more 
generally in open quantum systems [l9j . 

Recently, the description of open quantum systems 
by stochastic methods has received much attention (l9r - 



l2l| . In the Markovian limit, several methods have been 
proposed to treat fluctuation and dissipation starting 
from a quantum master equation on the system density 
[19 . |2(1 |22H28|| . These methods have been extended also 
to treat non-Markovian effects like in Quantum State Dif- 
fusion (QSD) [UrSi] or quantum Monte-Carlo (QMC) 
methods [331 ] . Several groups have shown that these 
effects could eventually be simulated using Feynman- 
Vernon influence functional [U, |34[ or directly stochastic 
master equations 35, 36]. 



In this work, we apply the stochastic formulation pro- 
posed in ref. [36] to the case of quadratic potentials 
coupled to a heat-bath, the so-called Caldeira-Leggett 
model 37]. The case of inverted potential is the first 
step towards realistic situations like fusion or fission. The 
aim of the present work is threefold. First, to introduce 
the new QMC method and apply it for potentials with 
barriers similar to those appearing in fusion/fission pro- 
cesses. Second, we show that the exact quantum Monte- 
Carlo method can be rather accurate to treat the dis- 
sipative dynamics of a quantum system. Last, we also 
present a comparison of this theory with other meth- 
ods based on projection, namely Nakajima-Zwangig (NZ) 
and Time-ConvolutionLess (TCL) [38M41I ] . which are ac- 
tually widely used to treat non-Markovian effects. Doing 
so, we show that only TCL up to at least fourth orders 
in the coupling constant can provide a competitive the- 
ory. The paper is organized as follows. In section [TTf the 
ingredients and properties of the quantum Monte-Carlo 
approach are discussed and the link with functional inte- 
gral is precised. In section [TTTl the method is first illus- 
trated to the case of parabolic potential. Then, the pass- 
ing probabilities are estimated for the inverted parabola 
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II. QUANTUM MONTE-CARLO METHOD 

Our starting point is a system (S) interacting with a 
surrounding environment (E). We assume here that the 
total system (S+E) is described by the Hamiltonian 

H = H S + H E + V. (1) 

Hs (resp. H E ) acts on the system (resp. env.) only 
while V induces a coupling between the two sub-systems. 
Starting from an initial total density -D(O), the dynamical 
evolution is given by the Liouville von- Neumann equation 
on the density: 

In many physical situations, the total number of degrees 
of freedom to follow in time prevents from solving exactly 
this equation. One of the leitmotiv of Open Quantum 
System (OQS) theory is to find accurate approximations 
for the system evolution without following explicitly irrel- 
evant degrees of freedom associated to the environment 
and therefore reduces the complexity of the initial prob- 
lem. Conventional strategy to treat dissipation and fluc- 
tuation in an open quantum system is to reduce the in- 
formation to the system density only ps(t) — TrE(D(i)) 
while accounting approximately for environment effect. 
Here, we use a different strategy, the dynamics of the 
total system is first replaced by a set of stochastic evo- 
lutions where the total density remains separable along 
each path, i.e. D = ps(t) <8> PE(t)- Then, the stochastic 
evolution of the environment is projected onto the rele- 
vant degrees of freedom to obtain a closed equation for 
the system density. It is shown that the new approach 
provides a proper treatment of dissipation and fluctua- 
tion for a system coupled to a surrounding heat-bath. 

A. Quantum Monte-Carlo formulation of Open 
Quantum Systems 

Recently, new stochastic formulations [H], [H, 52, 53 
have been developed to study the system+environment 
problem that avoid evaluation of non-local memory ker- 
nels although non-Markovian effects are accounted for 
exactly (see also [H, HI, E3, BE!)- One example of 
such a theory based on quantum Monte-Carlo method is 
presented here. 

Hereafter, it is assumed that the coupling is separable: 
V = Q ® B where Q and B act on the system and en- 
vironment respectively. For simplicity, initial separable 
density is considered, i.e. D(0) = ps{0) ® Pb(0). We 
want to replace the evolution of the total density (Eq. 

by an ensemble of stochastic evolutions of both the 
system and environment such that: 

dp s = § l [H s ,ps}+d^sQps + dX s psQ 

(3) 

d PE = §, [ h e, Pe] + di E Bp E + d\ E pEB 



where d^s/E an d dXs/E denote Markovian Gaussian 
stochastic variables with zero means and where we use 
Ito convention of stochastic calculus [48]. In the follow- 
ing, we assume in addition that 



dtsd\ E = dX s d^ E = 0, (4) 

where the overline denotes the stochastic average. Along 
each path, the total density remains separable, i.e. 
D(t) = ps{t) ® PE{t)- Starting from such a density, at 
time t + dt, the average evolution deduced from Eq. ([3]) 
is given by 



dD{t) = -[h s + h El D{t)\ 
in 

+dJsd^E (Q ® B)D(t) + dX s dX E D(t){Q ® B). 
Therefore, under the condition 



disd^E = tt, dX s dX E — ttt , (5) 
in in 

the average evolution over the separable densities that 
evolve according to Eq. ^ identify with the exact Li- 
ouville von Neumann equation of motion (Eq. (0)). The 
possibility to use simple Gaussian noises to incorporate 
the environment effect might appear surprising. Indeed, 
noises used in standard approaches for Open Quantum 
Systems generally reflect properties of the environment. 
It should be however kept in mind that such environment 
dependent noises appear once the environment dynam- 
ics has been projected out on the system density evo- 
lution. Anticipating the discussion of section III Bt once 
such a projection has been made, the Gaussian noises in- 
troduced here, transform into new random variables that 
explicitly depend on the environment properties. 

The discussion above for one time step can then be 
iterated to show that the exact dynamics of a sys- 
tem+environment could be replaced by an average over 
an ensemble of stochastic evolutions of separable densi- 
ties (34T4361 ]. To be really useful, mainly two difficulties 
should be overcome (i) in general, the environment cor- 
responds to a large number of degrees of freedom that 
could not be followed explicitly in time, (ii) the numeri- 
cal implementation of such a theory is possible only if the 
statistical errors do not growth too fast during the time- 
evolution. This statistical errors are directly connected 
to the number of trajectories necessary to accurately de- 
scribe the physical process. Difficulty (i) is solved in the 
next section by projecting out the effect of the environ- 
ment on the system leading to a closed equation for the 
system density only. Let us first concentrate on statis- 
tical errors. At any time, a measure of the statistical 
fluctuation around the average trajectory is given by 



A s4Qt = Tr{ (Z?t(i) - Dt(t)) (D(t) - D(t)) } 

= Tr{z?t£)(t)} -Tr{l3(t) 2 }. (6) 
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Starting from Eq. ([3]) , the evolution of A stat over a small 
time step reads 

dAstat = 2 -^{WT s +Wh}, (7) 

where <Q 2 ) g = Tr s (Q 2 p s (t)) and (B 2 ) E = 
TrE(B 2 j OE(t)). Statistical errors associated with Eq. ([3]) 
have been estimated numerically and turn out to grow 
very fast in time [46| . As a consequence, the stochastic 
process in the present form is useless to simulate physical 
situations and methods to reduce statistical errors should 
be used. 

To do so, it is worth to note that the stochastic equa- 
tion of motion is not unique. Indeed, any stochastic pro- 
cess of the form: 

' dp s = ft[H s + QA E , Ps ] 

+d£ s {Q - A s )ps + dXspsiQ - As) 

, (8) 

dp E = §[H E + BA s ,p E ] 

+d£ E {B - A E )p E + d\ EPE (B - A E ) 

where As(t) and A E (t) are time-dependent parameters 
leads to the same average evolution. These stochastic 
equations also provide a reformulation of the initial sys- 
tem+environment problem. Indeed, we have 



dps ® Pe + Ps® dp E = — [H s + QA E ,ps <g> p E ] 
in 

+% We + BA s ,p s ®Pe] 



dps® dp E = —[(Q-A s )®(B-A E ),ps®p E }. 
in 

Therefore, terms appearing in the deterministic part are 
exactly compensated by equivalent terms coming from 
the average over the noise. Accordingly, the evolution of 
the average density identifies with the exact equation of 
motion (j2|. 

Up to now, the flexibility has been essentially exploited 
by using [H, H|| 

A E (t) = (B(t)) E , A s (t) = {Q(t)) s . (9) 

This choice is justified by the fact that it directly appears 
when the Ehrenfest theorem is applied to separable to- 
tal density for system or environment observables. By 
modifying the stochastic evolution, part of the coupling 



is already contained in the deterministic evolution. Ac- 
cordingly, we do expect that the amount of coupling to 
be treated by the noise is significantly reduced as well 
as the statistical errors. In the latter case, statistical 
fluctuations are given by: 



/tit r 2 n i 

dKtat = -y{((Q 2 )s-(Q) s ) + ((B 2 )s~(B) s )}iO) 

and are always smaller than the original ones ([7]). As 
shown numerically in ref. (35[, introduction (j9|) signifi- 
cantly reduces statistical fluctuations and opens new per- 
spectives for the application of the present framework. 

The modified stochastic theory has other advantages. 
For instance, the traces of densities are constant and 
remain equal to their initial values, i.e. dTr(ps/ E ) — 
0. This greatly simplifies expectation values of system 
and/or environment observables. Indeed, denoting by X 
a system operator, along a trajectory, we have 

(X) = Trv(XD(t)) = Tr(p E (t))Tr s {X Ps (t)). (11) 

For stochastic processes with varying trace of densi- 
ties, the observable evolution will contain terms com- 
ing from dTr(ps/ E ) and cross terms coming from 
dTr(ps/ E )dTr(Xps(t)). In the case considered here, we 
simply have 

d(X) = Tr E (p E (t))dTr s (Xps(t)). (12) 

The QMC theory with centered noise overcomes the dif- 
ficulty (ii) but does not help for (i) since the environment 
degrees of freedom should still be followed in time. In the 
next section, we show how irrelevant degrees of freedom 
can be projected out to obtain a closed stochastic master 
equation for the system only. 



B. Reduced system density evolution and link with 
influence-functional theory 

The stochastic formulation suffers a priori from the 
same difficulty as the total dynamics: the environment 
is in general rather complex and has a large number of 
degrees of freedom which can hardly be followed in time. 
In Eq. |(5} , the influence of the environment on the sys- 
tem only enters through {B{t)) E . Therefore, instead of 
following the full environment density evolution, one can 
concentrate on this observable only. As shown in ref. 
|43| . the second equation in Eqs. © can be integrated 
in time to give: 



{B(t)) E = Trv{B I {t-t )p E {t ))-\ f D{t,s)(Q{s)) s ds- f D(t, s)du E {s) + f D 1 {t,s)dv E (s), (13) 

n Jo Jo Jo 
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where B 1 denotes the operator B written in the interaction picture while D and D\ are the memory function given 
by 

D(t,s) = i([B(t),B(t-s)]) E , and £>i(t, s) = ({B(t), B(t - s)}+) E - 2(B(t)) E {B(t - a)) E . (14) 

A new set of stochastic variables dvs/ E and dus/E have been introduced through d£s/E — dv$/E — idus/E an d 
dXs/E = dv s/E + idu s/E , and verify 



dusduE — dvsdvE = rr, dugdvE — dvsduE — 0. (15) 

Reporting the evolution of (B(t)) E into the evolution of p$, a closed stochastic equation of motion for the system 
density is obtained: 

dt 

d Ps = ^[H s ,Ps] + dt[Q,p s ] / dsD(t-s)(Q(s)) s + d^t)[Q,p s ]+dr}(t){Q-(Q) s ,p s } (16) 

I'll Jo 

with 

d£(t) = dt ( Di(t- s)dv E {s) -dt ( D(t - s)du E (s) - idv s (t), drj{t) = du s {t). (17) 
Jo Jo 



By integrating out the evolution of the environment, a 
new stochastic term is found that depends not only on 
the noise at time t but also on its full history through 
the time integrals. Using second moments given by Eqs. 
(fl~5j) leads to: 

dr)(t)dr)(t') = 0, 

dW)dW) = ~Q(t-l/)D(t-l/), 
dWdW) = ~D x {\t-H\), 

where <d(t — t 1 ) = 1 if t > t' and elsewhere. Inter- 
estingly enough, the stochastic equation given by (flB)) 
identifies with the stochastic master equation obtained 
in ref. [211 ] using a completely different method based on 
the Feynman- Vernon influence functional theory [49(. It 
should however be kept in mind that different strategies 
to design the stochastic equation (see discussion in sec- 
tion III A|) would have given a different stochastic master 
equation. 

Despite the apparent complexity of Eq. (fT6|) , the QMC 
approach has been recently applied with success to the 
spin-boson model coupled to a heat bath of oscillators 
[431 ] . In particular, the introduction of Q seems to cure 
the numerical difficulties that have been encountered in 
this model [46| . By projecting the environment effect 
onto the system density evolution we do not need any- 
more to follow the environment density and we expect 
that a rather limited number of trajectories will be suffi- 
cient to accurately simulate the onset of dissipation and 
fluctuation in an open quantum system. Eq. (fTB)l is 
the equation that is solved in practice. It should be 
noted that the present stochastic process differs signif- 
icantly from conventional approaches. Indeed, accord- 
ing to the noise properties, system densities are non- 



hermitian along a stochastic path. As a consequence, ex- 
pectation value of observables are complex loosing their 
physical meaning before averaging over the stochastic 
path. Nevertheless, we illustrate in the following that 
the new theory can be a very powerful tool. 

III. APPLICATION 

The Caldeira-Leggett (CL) model [37| corresponds to 
a single harmonic oscillator coupled to an environment of 
harmonic oscillators initially at thermal equilibrium, i.e. 

H s = H c + ^+e^Mu 2 Q 2 , (18) 
H E = E + \m n ^l) (19) 

and B = -J2n K nX n HH- Here, H c = Q 2 J2„ 2 m„"^ is 
the counter-term that insures that the physical frequency 
is luq. In the following, e is cither +1 (harmonic case) or 
— 1 (inverted parabola case). Such a model can be solved 
exactly. 

As shown in ref. |43| . the two functions D and Di 
defined by Eqs. (fl"4")l and estimated along the stochastic 
trajectories identify with the standard times correlation 
functions: 

r + oo 

D(t) = 2h dojJ(u})sm(ojT), (20) 
Jo 

r + oo 

Di(t) = 2h duj J(u) coth(Huj/2kBT) cos(cjt)(21) 
Jo 

K 2 

where J(cj) = - — - — 8 (to — w„) denotes the spectral 
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density characteristic of the environment [19|, |50j ■ In the 
following, a Drude spectral density [43[ 



ft 2 



ft 2 ' 

is considered where M is the nucleon mass. 



(22) 



A. Quantum Monte-Carlo method for parabolic 
potentials 

As a first illustration, the harmonic case (e = +1) is 
considered. This case has been already studied in ref. 
[oil ] using the stochastic method proposed in ref. [2l[ . In 
the CL model, starting from an initial Gaussian density, 
the system density remains gaussian along the stochastic 
path. Therefore, the stochastic evolution of the system 
density reduces to the first and second moments evolution 
of (P) and (Q) given by [H: 



d(Q) = ^dt + 2du s a QQ 

d (P) = -Mu>l{Q)dt - dt (B) + 2du s a PQ - Mv s 



d<JQQ 

da pp 



2 W a PQ 
—ZMw^dtapQ 



(23) 



da 



PQ 



a PP - MwlaQQdt 



These equations illustrate the differences between the 
new exact reformulation and standard methods to treat 
dissipation. Generally, the noise enters into the evolu- 
tion of (P) only and affects directly the second moment. 
Here, we see that second moments identify with the un- 
perturbed ones while the random forces enter in both (Q) 
and (P). In addition, the noise is complex, which im- 
plies that observables make excursions into the complex 
plane. This stems from the specific noise used to design 
the exact formulation that leads to non-Hermitian den- 
sities along paths. Part of the conceptual difficulties in 
understanding the physical meaning of observable evolu- 
tions can be overcome by noting that if ps(t) belongs to 
the set of trajectories, by symmetry p s (t) will also be- 
longs to the set. By grouping these two trajectories to 
estimate observables, real quantities are deduced. 

The exact evolution is obtained by averaging over dif- 
ferent trajectories. For second moments, this leads to 



Eqq = (Q 2 ) - (Q) = a QQ + (Q) 2 - (Q) 



(24) 



where (X) denotes the statistical average of quantum 
expectation values (X). It is a particularity of the 
CL model that total fluctuation is recovered by simply 
adding up quantum and statistical fluctuations. 

An example of Sgg(t) evolution obtained using Eq. 
(IM)l (red filled circles) is compared to the exact result 
(solid line) in figure [TJ As an indication, the evolution of 




FIG. 1: (Color online) Evolution of Eqq (filled circles) ob- 
tained by averaging over 10 5 trajectories. This evolution is 
compared to the exact result (solid line) and to the quan- 
tum fluctuation oqq evolution (dotted line). Parameters of 
the simulation are UbT = fkoo, Ml = 5huo, r\ — 0,5frwo and 
hw = lAMeV. The factor 7 defined as -y 2 = fi/(2Mw ) 
equals here 7 = 1.216 fin - . 



quantum fluctuation aQQ, which is identical for all tra- 
jectories, is also displayed (dotted line). Note that, Eq. 
(I16p is already exact for the evolution of first moments 
(P) and (Q), even if the noise is omitted. However, it 
completely fails to account for fluctuation. While the 
quantum evolution does not present any damping, the 
average evolution closely follows the exact solution. The 
harmonic oscillations in aQQ are due to the fact that the 
width of the initial density differs from the width of the 
coherent state associated to the considered oscillator, i.e. 
ctqq(O) =/= h/(2AluJo). This is at variance with the simu- 
lation made in ref. [5l| . 

The accuracy of the quantum Monte-Carlo theory has 
been systematically investigated for various temperatures 
and coupling strengths. In all cases, averaged evolutions 
could almost not be distinguished from the exact evo- 
lution. This is illustrated in figure [5] where Sqq (left), 
Epp (middle) and Epq (right) are displayed as a func- 
tion of time and compared to exact solutions for various 
temperatures. Figure [21 clearly shows that the stochas- 
tic method properly includes all non-Markovian effects. 
In particular at low temperature, typically kpT < Hloq, 
and medium coupling constant 77, large memory effect is 
expected. 



B. Application of NZ and TCL 

Conjointly to the benchmark of quantum Monte-Carlo 
approaches, we also tested projection method either 
based on the Nakajima-Zwa nzig |19l . l38l . 139 . |52| or Time 
ConvolutionLess [H, E3> EL Hlf formalisms. Both theo- 
ries provide a priori exact re-formulations of the initial 
problem and lead to a closed master equation for the sys- 
tem density. However, they differ completely in the strat- 
egy and equation of motion used to incorporate memory 
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FIG. 2: (Color online) Evolution of E PP (left), Eqq (middle) 
and Epg (right) obtained with 10 s trajectories are displayed 
with red filled circles as a function of time and systematically 
compared with the exact evolution (solid line). ksT = 5huo, 
fkoo and 0, 1/kJo are respectively shown from top to bottom. 
In all cases, rj — 0, 5&Jo, Ml = 5fio>o and huuo = 14MeV. 



effects. In the NZ case, the evolution of the system den- 
sity at time t depends on its full history (i.e. on ps(s) 
for all s < t). In the TCL case, the master equation is 
local in time and non-Markovian effects are treated in 
time-dependent transport coefficients. To illustrate the 
differences between our new QMC method and TCL, we 
remind below the corresponding local master equation, 
for V = Q ® B. 



t d , s 
K Jt p s {t) 



i[H Sl p s (t)}--A(t)[Q,{Q,p s (t)}} 
2iX(t)[Q, {P, Ps (t)}] - %^[Q, [Q, ps(t)}} 



■2?qV[Q,[p,p s (t)}]. 



(25) 



A(t), X(t), DpQ{t) and Dpp{t) are time-dependent 
transport coefficients that contain memory effects. Simi- 
larly to the QMC case, the solution of the master Eq. 
(I25p is equivalent to follow first and second moments 



given by: 

r <wi 

dt 

d(P) 
dt 



dE 



pp 



dt 

d^QQ 



dt 



dE 



PQ 



dt 



M 

= - Mo£(t) (Q)-2X(t)(P) 

= - 2Mu3 2 p {t) E PQ -4A(i)E PP + 2D PP (t) 

- o Sp Q 
M 

= - Mu 2 p (t) E QQ -2A(i)Sp Q 



M 



2D PQ {t) 



with Up(t) = ujq + A(t). In practice, the exact NZ or 
TCL theory cannot be exactly solved and an expansion in 
powers of the coupling constant is made. In the following, 
NZ2 (or TCL2) will refer to the expansion up to second 
order while NZ4 (or TCL4) will refer to the expansion is 
made up to fourth order. By neglecting higher orders in 
the coupling in NZ2 (resp. TCL2) or NZ4 (resp. TCL4), 
both theory are not exact anymore. In the following, the 
efficiency of each method is systematically discussed. 

In Figure |31 the evolution of Epp for different cut-off 
frequencies Ml and coupling strengths rj are compared 
to the exact evolution (solid line). Explicit forms of the 
equation of motion in the NZ and TCL case can be re- 
spectively found in ref . [H, [H, H, and [H, HJ HE . 
Several important remarks could be drawn from this com- 
parison: (i) In all cases, when the coupling strength is 
considered up to second order, NZ2 (open triangles) pro- 
vides a better approximation than TCL2 (open squares) . 
This might indeed be expected since NZ2 and TCL2 are 
respectively equivalent to the Born and Redfield master 
equation and the former contains a priori less approxima- 
tions than the latter, (ii) While the TCL4 (filled squares) 
leads to a clear improvement compared to the TCL2, NZ4 
(filled triangles) is in general worse than NZ2. This is a 
known difficulty of NZ approach and was one of the mo- 
tivation for the introduction of TCL method (see discus- 
sion in ref. H3, EH). This stems from the fact that the 
order in perturbation in NZ cannot be identified. For in- 
stance, NZ2 (resp. NZ4) contains orders in coupling con- 
stant greater than 2 (resp. 4). As a result, the NZ theory 
does not lead to better results when the " apparent" order 
in the coupling increases. The TCL method essentially 
cures this pathology and precise orders in the coupling 
can be selected, (iii) Rather large deviations between 
the exact and TCL2 are observed for different cut-off fre- 
quencies and coupling strengths. This issue is important 
since several theories have been recently developed along 
the line of TCL2 to include memory effects in fusion and 
fission reactions [3, EEI, E3] • Note that, the accuracy of 
TCL2 depends on different parameters used in the spec- 
tral density, in particular of the parameter Ml. Here, 
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we have used a value of the cut-off frequency between 10 
and 20 MeV, which gives realistic dissipation and fluc- 
tuation for the fusion or fission mechanism [13, EH, EH- 
Our study clearly points out that a proper treatment of 
memory requires to include higher order effects, (iv) Fi- 
nally, in all cases, TCL4 could not be distinguished from 
the exact result. As we will see, the efficiency of TCL4 
is similar for the inverted parabola. 




FIG. 3: (Color online) Evolution of Epp for different approx- 
imations: NZ2 (open triangles), NZ4 (filled triangles), TCL2 
(open squares) and FCL4 (filled squares). The exact evo- 
lution is displayed with solid line. In all cases, huoo = 14 
MeV, and ksT = hwo are used. The left side, corre- 
sponds to different cut-off frequencies: HQ — 20fiu>o (top) and 
hQ. — 5hwo (bottom). In both cases, r\ — 0, bhuo- In the right 
side, hfl = lOfioJo and different coupling strengths are used: 
n = huio (top) and r\ = 0.1/iwo (bottom). 

Since NZ method is not competitive, only the quan- 
tum Monte-Carlo and TCL methods are considered in 
the following application. 

C. Quantum Monte-Carlo method applied to 
inverted oscillators 

Several approaches have been recently developed to de- 
scribe fusion and fission reactions @, EH, EE E3, El, HH] • 
In these mechanisms, few collective degrees of freedom 
couple to a sea of internal excitations while passing an 
inverted barrier. At very low energy, both quantum and 
non-Markovian effects are expected to play a significant 
role. Most of the theory currently used start from quan- 
tum master equations deduced from TCL2. The quan- 
tum Monte-Carlo method offers a practical alternative 
which has similarities with path integrals theory. Path 
integrals are known to provide a possible framework to 
include dissipation while passing barriers (see for instance 
[54j). However, due to their complexity, only few appli- 
cations have been made so far [2, 55]. We compare here 
the different approaches for inverted potential (e = — 1). 



1. Initial conditions, trajectories and mean evolution 

Initially, we consider a Gaussian density with quantum 
width ctqq(0) = 0.16 fm 2 and ctpq(0) = MeV.fm/c 
and positioned on one side of the potential (here taken 
arbitrarily at (Q{0)) — Qo > while the barrier height is 
is located at fm and is by convention taken as Vb = 
MeV). The initial kinetic energy, denoted E K (0) is set by 
boosting the density with an initial momentum (P(0)) = 
Po < 0. 

Contrary to the classical theory of Brownian motion, 
the notion of trajectories is not so easy to tackle in the 
present Monte-Carlo framework. First, observables are 
complex. As mentioned in section MI Al this difficulty 
can be overcomed by grouping trajectories by pairs which 
is equivalent to replace expectation of observables by 
their real parts. Second, it should be kept in mind that 
the present theory is a purely quantum theory where den- 
sities associated to wave-packets are evolved. Therefore, 
each trajectory should be interpreted in the statistical 
sense of quantum mechanics and contains many classical 
paths. Nevertheless, to visualize the trajectory we define 
the following energies: 

where Q(t) and P{t) denote the real part of (Q(t)) and 
(P(t)) along the trajectory. An illustration of two trajec- 
tories, one passing the barrier and one reflected is shown 
in figure |U As illustrated in the following, it is conve- 
nient to group trajectories according to the quantity AE 
defined by 

AE = £7(0) - V B (27) 

which is nothing but the difference between total initial 
energy and barrier high. Both trajectories shown in fig- 
ure 0] correspond to AE = MeV. 

It is tempting to group trajectories into those passing 
the barrier and those reflected by the potential to get 
information on the passing probability or passing time, 
however, it should be kept in mind that the present the- 
ory is fully quantal. Since each trajectories are associ- 
ated with densities with quantum widths, both trajecto- 
ries presented in Fig. @] contribute to the transmission 
probability. 

The accuracy of different methods is illustrated in fig- 
ure [5] where evolutions of (Q), (P), Eqq and £pp are 
shown as a function of time. Values of parameters re- 
tained for this figure are typical values generally taken in 
the nuclear context [IB] . In all cases, including TCL2, sec- 
ond moments are well reproduced. However, only TCL4 
and the stochastic simulation provides a correct descrip- 
tion of first moments. Calculations are shown here for 
AE = MeV. TCL2 provides a better and better ap- 
proximation when AE increases while the disagreement 
increases below the barrier. This will be further illus- 
trated below. 



D. Transmission probability 
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FIG. 4: (Color online) Evolution of E(t) as a function of Q(t) 
for two trajectories (dark and bright lines) with A_E = MeV, 
k B T = lMeV and rj = 0.003MeV. The black arrow indicates 
the initial position of the trajectories while the potential is 
also shown with bold line as a reference. 

Different coupling parameters, cut-off frequencies and 
temperatures have been investigated showing that both 
TCL4 and quantum Monte-Carlo are very accurate the- 
ories leading always to results on top of the exact ones. 
It should be noted that the number of trajectories used 
in the stochastic approach to get small statistical errors 
is rather small (around 10 5 ) which is quite encouraging 
for future applications. 
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FIG. 5: (Color online) Evolution of (Q), (P), £qq and Epp as 
a function of time obtained with quantum Monte-Carlo (filled 
circles), TCL2 (open squares) and TCL4 (filled squares). The 
exact evolution is also displayed with a solid line. Param- 
eters of the simulations are Vb = 4MeV, fiWQ = lMeV, 
rj = 0.03MeV, Ml = Ihhu), k B T = IMeV, AE = MeV 
and 7 = 0.402 fm . 10 5 trajectories have been used for the 
quantum Monte-Carlo case. 



The accuracy of the method used to incorporate non- 
Markovian effects directly affects the predicting power of 
the theory. This aspect is illustrated here with the pass- 
ing probabilities. Such a probability is a crucial ingredi- 
ent in particular for models dedicated to the formation of 
very heavy elements @, 0, EE EE HE HE H3] and should 
be precisely estimated. 

The asymptotic passing probability is usually defined 
as: 

P(+oo) = lim ierfc f ] . (28) 

where q(t) and cr qq {t) denote the expectation value and 
second moment of Q deduced from the considered the- 
ory. In the quantum Monte-Carlo case, these quantities 
identify with (Q(t)} and £qq(£) respectively. To quan- 
tify the precision of each theory, we have systematically 
investigated the difference between the estimated pass- 
ing probability and the exact one using the parameter 
AP/P: 

AP _ P(+oo) - Pex(+OQ) 
P ~ P cx (+oo) 

where P(+oo) and P ex (+oo) denote the results of the 
specific calculation considered and the exact one respec- 
tively. Figure [S] presents the evolution of AP/P as a func- 
tion of AE obtained for different coupling strengths and 
temperatures for the quantum Monte-Carlo (filled cir- 
cles), TCL2 (open squares), TCL4 (filled squares) cases. 
In this figure, the Markovian approximation is also shown 
by open crosses. 

Once again, the TCL4 and the quantum Monte Carlo 
methods are in perfect agreement with the exact solu- 
tion for any input parameters. Small differences some- 
times observed between the stochastic approach and the 
exact value come from the limited number of trajectories 
used to obtain figure |S] Well above the barrier (here at 
least two times), TCL2 converges towards the exact case. 
However, at low AE, it turns out to be a rather poor ap- 
proximation. The difference seen in the TCL2 case can 
directly be traced back to the discrepancy already ob- 
served in figure [5] and further confirms the difficulty of 
treating non-Markovian effects below the barrier. We can 
see that at lowest energy considered here, the error could 
be as large as 20% in the weak coupling case and more 
than 100% in the strong coupling limit. It is worth finally 
to mention that below the barrier, the Markovian limit 
gives an even better result that the TCL2 case. 

IV. SUMMARY AND DISCUSSION 

In this work, the quantum Monte-Carlo approach re- 
cently proposed in ref. [43] to incorporate exactly non- 
Markovian effects is introduce and applied to the case of 
harmonic potentials coupled to a heat-bath. 
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FIG. 6: (Color online) Evolution of AP/P as a function of 
AE calculated using quantum Monte-Carlo (filled circles), 
TCL2 (open squares), TCL4 (filled squares) and Markovian 
approximation (open crosses) for different coupling constant: 
t] = 0.003MeU (left) and rj = O.OSMeV (right). In both 
cases, T = 5hw, hui and 0, lJuo are shown from top to bottom. 



Other methods have also been benchmarked. The TCL2 
method, which is now widely used in nuclear physics 
to estimate passing probabilities, turns out to deviate 
significantly from the expected result especially below 
the barrier and even in the weak coupling regime. To 
properly treat the dynamics of barrier transmission, 
higher orders in the coupling strength should be incorpo- 
rated. TCL4 gives very good agreement with the exact 
evolution in all cases considered here. The conclusion of 
our present work is that both quantum Monte-Carlo ap- 
proach and TCL4 could be considered as good candidates 
to include memory effects for situations of interest in 
nuclear physics. Henceforth, the TCL2 method which is 
widely used nowadays should be replaced by TCL4. The 
possibility to use stochastic formulation that are exact 
in average open new perspectives to describe a system 
coupled to a complex environment. The application to 
harmonic potential gives interesting insight into such a 
theory. Application to more general potential turns out 
to be less straightforward with the appearance of spikes 
which have been already observed in several formalism 
where non linear stochastic equations appear [28[. To 
make these theories more versatile, new methods like 
the one proposed recently in ref. |58j | could be used. 
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